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We study semiclassical dynamics of anisotropic Heisenberg models in two and three dimensions. 
Such models describe lattice spin systems and hard core bosons in optical lattices. We solve numer- 
ically Landau-Lifshitz type equations on a lattice and show that in the phase diagram of magnetiza- 
tion and interaction anisotropy, one can identify several distinct regimes of dynamics. These regions 
can be distinguished based on the character of one dimensional solitonic excitations, and stability 
of such solitons to transverse modulation. Small amplitude and long wavelength perturbations can 
be analyzed analytically using mapping of non-linear hydrodynamic equations to KdV type equa- 
tions. Numerically we find that properties of solitons and dynamics in general remain similar to 
our analytical results even for large amplitude and short distance inhomogeneities, which allows 
us to obtain a universal dynamical phase diagram. As a concrete example we study dynamical 
evolution of the system starting from a state with magnetization step and show that formation of 
oscillatory regions and their stability to transverse modulation can be understood from the proper- 
ties of solitons. In regimes unstable to transverse modulation we observe formation of lump type 
solutions with modulation in all directions. We discuss implications of our results for experiments 
with ultracold atoms. 



Motivation. Understanding nonequilbrium quantum 
dynamics of many-body systems is an important prob- 
lem in many areas of physics. The most fundamental 
challenge facing this field is identifying emergent col- 
lective phenomena, which should lead to universality 
and common properties even in systems which do not 
have identical microscopic Hamiltonians. We know that 
in equilibrium many-body systems typically fall within 
certain universality classes. Basic examples are states 
with spontaneously broken symmetries such as BEC of 
bosons and superfluid paired phases of fermions[l, 2], 
Fermi liquid states of electrons and interacting ultracold 
fermions[3, 4]. We have numerous examples of emer- 
gent universal properties of classical systems driven out 
of equilibrium (see ref. [5] for a review). However very 
few examples of collective behavior of coherent quantum 
evolution of many-body systems, which can be classified 
as exhibiting emergent universality, are known (several 
recently studied examples can be found in refs [6-12]) 
The main result of this paper is demonstration of the 
universality of semiclassical dynamics of lattice spin sys- 
tems and strongly interacting bosons in two and three 
dimensions, summarized in the phase diagram in Fig. 1. 

In addition to fundamental conceptual importance, 
our study has strong experimental motivation. Rapid 
progress of experiments with ultracold atoms in optical 
lattices makes it possible to create well controlled real- 
izations of systems that we discuss [13, 14]. Tunability of 
such systems, nearly perfect isolation from the environ- 
ment, and a rich toolbox of experimental probes, includ- 
ing the possibility of single site resolution [15-18], makes 
them excellent candidates for exploring nonequilibrium 
quantum dynamics (see [10, 19] and references therein). 
Several intriguing nonequilibrium phenomena have been 
demonstrated recently in ensembles of ultracold atoms, 
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FIG. 1: Dynamical phase diagram of the anisotropic Heisen- 
berg model (1) for initial states with finite magnetization in 
the XY plane. Vertical axis labels magnetization for spin 
models and density for hard-core bosons on a lattice. Re- 
gions I, II correspond to easy plane anisotropy of the Hamil- 
tonian, linearized collective modes are of the hyperbolic type. 
Regions III correspond to easy axis anisotropy, linearized col- 
lective modes are of the elliptic type. Special lines J z /J± = 1 
and m z — ±1 all have linearized collective modes of parabolic 
type. Region 1+ has particle solitons unstable to transverse 
modulation. Region I_ has hole solitons unstable to trans- 
verse modulation. Region 11+ has hole solitons stable to 
transverse modulation. Region I_ has particle solitons sta- 
ble to transverse modulation. 



including collapse and revival of coherence [20], absence 
of relaxation in nearly integrable systems [6], exponen- 
tial slowdown of relaxation in systems with strong mis- 
match of excitation energies[21], anomalous diffusionfll], 
prethermalization[12], many-body Landau-Zener transi- 
tions [22] , light-cone spreading of correlations in a many- 
body system [23]. 

Model. In this paper we study non-equilibrium dynam- 
ics of lattice spin models and strongly correlated bosons 
in optical lattices. Our starting point is the anisotropic 
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Hcisenbcrg model 

H AH = -Jx E + -J'Y, W 

where cr a arc Pauli matrices. Anisotropic Hcisenbcrg 
models with tunable interactions can be realized with 
two component Bose mixtures in spin dependent opti- 
cal lattices [24-27]. Hamiltonian (1) also describes spin- 
less bosons in the regime of infinitely strong on-site 
rcpulsion[28]. In this case states | 1) and | t) correspond 
to states with zero and one boson per site respectively, 
J± is the tunneling strength, and J z is the strength of 
the nearest neighbor repulsion. For atoms with contact 
interaction and confined to the lowest Bloch band J z = 0. 
Non-local interactions can be present for atoms in higher 
Bloch bands [29] and polar molecules [30]. 

Semiclassical equations of motion can be understood 
either as the lattice version of Landau-Lifshitz equa- 
tions or as dynamics in the manifold of variational states 
!*(*)> = n> in J e"W2| ^ +C0S | e W2| ^ whcrc 
all 0i and tpi are independent functions of time ( see sup- 
plementary material and ref. [31]). It is convenient to 
introduce parametrization 

sin | = (i^) 1 / 2 , cosf = CA^) 1 ' 2 , e** = 
,2^2^1/2 ; then semiclassical dynamical equations arc 
given by 

±i = 4J±z l (y i+ i + yi_x) ~ <±JzVi (zi + i+Zi-i) 
Vi = -4 J±Zi (x i+1 + 4 J z Xi (Zi+i+Zi-i) 

Zi = 4 Jj_ {y l x l+1 - x,y i+ i + i/.i,^ - 5CiJ/»-i) 

(2) 

Note that variables (xi,yi,Zi) reside on a sphere. Speci- 
fying xf + yf + zf =1 in the beginning of the evolution 
will preserve this condition at all times. 

We consider dynamics in a state with finite magneti- 
zation in the XY plane (supcrfluid phase for hard core 
bosons), where close to equilibrium dynamics is deter- 
mined by the Bogoliubov (Goldstone) mode. Consider- 
able difference between the nonlinear hydrodynamics of 
hard core bosons on a lattice and the more familiar GP 
model was observed in the numerical study of solitary 
waves in ref. [32] (for a discussion of solitons in systems 
of ultracold atoms in the regime where GP equation ap- 
plies see refs. [33-44]). Additional evidence for the spe- 
cial character of solitons in systems of strongly corre- 
lated lattice bosons came from the analytical study in 
Ref. [31], which extended linear hydrodynamics to in- 
clude the first non-linear term and the first dispersion 
correction to the long wavelength expansion and showed 
that solitons were of the KdV type. In this paper we 
derive a full phase diagram of semiclassical dynamics of 
(1) for initial states with finite magnetization in the XY 
plane. Our approach relies on combining analytical re- 
sults describing small density perturbations [31] with nu- 
merical studies of equations (2) describing dynamics of 
large amplitude modulations. 



Phase diagram. Hyperbolic, elliptic, and parabolic 
regimes. The main results of our analysis are summa- 
rized in fig. 1. Firstly there is fundamental difference 
between elliptic and hyperbolic regions in the character 
of linearized equations of motion. The hyperbolic regime 
(hyperbolic regions I and II in fig. 1) corresponds to the 
easy plane anisotropy of the Hamiltonian. If we start 
with initial state that has finite magnetization in the XY 
plane, we find collective modes with linear dispersion, 
u) oc \q\, which describe Bogoliubov (Goldstone) modes. 
In the elliptic regime we have easy axis anisotropy: lowest 
energy configuration should have full polarization along 
the z-axis, m = ±1. However our initial state has mag- 
netic order in the XY plane. Thus we find unstable collec- 
tive modes with imaginary frequency Imtu oc i\q\. This 
corresponds to the dynamical instability in which small 
density fluctuations grow in time exponentially at short 
time scales. There are also several special lines on the 
phase diagram 1: the SO (3) symmetric Heisenberg model 
with J z = Jj_ and fully polarized regimes m = ±1 (fully 
occupied or empty regimes for spinless bosons). In all of 
these cases collective modes have quadratic dispersion, 
i.e. linearized equations of motion are of the parabolic 
type. 

Hyperbolic regime. Solitons. Hyperbolic regions I and 
II arc differentiated further according to the character of 
solitonic excitations. In region I + we find particlc-like 
one dimensional solitons, which arc unstable to modula- 
tion in the transverse directions. In region I_ we find 
hole-like one dimensional solitons unstable to transverse 
modulations. In regions II± we find particle and hole-like 
solitons stable to transverse modulation. Lines IV i and 
IV2 correspond to two different regimes of the mKdV 
equation. Both hole and particle-type soliton solutions 
can be found on the line IV2 . Both types of solitons are 
suppressed close to the line IVi. 

In the limit of small amplitude and long wavelength 
inhomogencitics non-linear hydrodynamic equations (2) 
can be mapped to KdV type equations and solitons can 
be studied analytically (see supplementary material and 
ref. [31]). This analysis can be used to obtain explicit ex- 
pressions for boundaries of different regimes in fig. 1. In 
the general case of large amplitude inhomogeneities one 
needs to solve lattice equations (2) numerically. Surpris- 
ingly we find that in even for large amplitude solitons, 
their properties, including stability to transverse fluc- 
tuations, arc consistent with small amplitude solitons. 
Hence the phase diagram in fig. 1 is generic. 

While we can not present all of our numerical results, 
we provide examples, which demonstrate properties of 
large amplitude solitons. Fig. 2 shows particle- and hole- 
type solitons in regions 1+ and 11+ respectively. These 
specific solutions were chosen arbitrarily from a plethora 
of solitons with different amplitudes and velocities, which 
can be found in the system for the same values of mi- 
croscopic parameters. In these specific examples only 
one dimensional variation of parameters were considered 
when solving lattice equations (2). The possibility of 



FIG. 2: Left figure: particle soliton in region I+. Right fig- 
ure: hole soliton in region 11+ . Red line shows magnetization 
(density) cos ft. Blue line shows superfluid velocity <pi+i — <{>i. 
Results obtained from numerical solutions of lattice equations 
(2). Transverse modulation was not included. 




FIG. 4: Time evolution of a soliton for J z /J± = in region 1+ 
in a 2d system. Axis are the same as in fig. 3 Different plots 
correspond to t J± = 0, 40, 80, 90. This soliton is unstable 
to small transverse modulation introduced explicitly in the 
initial state. 




FIG. 3: Time evolution of a soliton for J z /J± = in region 
11+ in a 2d system. Horizontal axis are spatial coordinates. 
Vertical axis is magnetization (density). Different plots cor- 
respond to tj±_ = 0, 10, 30, 50. This soliton is stable to 
to small transverse modulation introduced explicitly in the 
initial state. 



transverse modulation of one dimensional solitons in two 
dimensional systems is considered in figures 3, 4. These 
figures contrast stable and unstable regimes. In the sta- 
ble regime II + transverse modulation does not lead to 
any dramatic change of the soliton solution. Analo- 
gous dynamics takes place in the other stable regime, 
II_, except for a change from hole solitons to particle 
ones. In the unstable regimes 1+ (the same behavior 
is observed in I_) we observe that transverse modula- 
tion leads to the formation of two-dimensional "lump" 
solutions. For small amplitude modulations stability to 
transverse modulation can be studied analytically using 
Kadomtsev-Petviashvili equation [31, 45]. Our analysis 
of the stability to transverse modulation is not limited to 
solitons but applies to all inhomogencous states. 

Parabolic and elliptic regimes. A rapidly growing os- 
cillation zone characterizes elliptic instability for system 
(2) in region III on the phase diagram. It is shown in 
fig. 5. In supplementary material we show solution of 
(2) in the parabolic regime separating the elliptic and 
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FIG. 5: Time evolution of a small density perturbation in 
the elliptic region. Vertical axis is magnetization (density). 
Note rapidly developing oscillation zone typical to modulation 
instability. Different plots correspond to tj± = 0, 2, 4, 6. 



hyperbolic regions on the phase diagram 1. 

Decay of magnetization step. We now discuss the prob- 
lem of the decay of magnetization step in a system of 
type (1). Configurations of this type have been recently 
realized in 3d magnetic systems using tunable magnetic 
field gradient [46]. They can be created for strongly cor- 
related spinless bosons in 2d optical lattices using local 
addressability [17, 18]. Besides this experimental motiva- 
tion, relaxation of the magnetization (density) step is an 
important methodological problem. Dynamics starting 
from this state allows extended unattenuated propaga- 
tion of the magnetization (density) modulation, which 
should amplify the role of nonlinearities. This may be 
contrasted to e.g. dynamics starting with a ring type in- 
homogeneity, where already at the level of linear hydro- 
dynamics there is a decrease of the magnetization (den- 



FIG. 6: Decay of the density step in Region I+. Red lines de- 
scribe density, blue line corresponds to the phase difference, 
(a) Initial state t — 0. (b) After t J± = 25 of the evolution. 
Note separation of the right and left moving parts and for- 
mation of the oscillatory region on the left moving part, (c) 
Expanded view of the left front after tj± = 5000 showing 
formation of the leading edge of the oscillation zone, (d) Left 
front after t.J±_ = 15000. Note soliton pairing mediated by 
background waves with small amplitude. 



sity) modulation in time. In classical hydrodynamics de- 
cay of the density step is one of the canonical problems 
considered in Refs. [47-49]. 

Fig. 6 shows the main stages of the magnetization step 
decay in regions I± , II± , IV: separation of left- and right- 
moving parts, steepening of one of the moving edges and 
formation of the oscillatory front. While there is general 
agreement between small amplitude limit analyzed in ref. 
[31] and results of numerical analysis of lattice equations 
presented in fig. 6, there is one important difference. In 
the oscillatory region appearing from the decay of a large 
amplitude magnetization step we observe pairing of soli- 
tons, which can be understood as a result of interaction 
between solitons through the background of small ampli- 
tude waves. This is manifestation of the absence of exact 
intcgrability of lattice equations (2). 

Experimental considerations. Direct experimental 
tests of theoretical results presented in this paper can 
be obtained by analyzing dynamics starting from soliton 
configurations of ultracold atoms in optical lattices, such 
as shown in figs 2, 3, 4. These figures show that lat- 
tices of the order of tens of sites should be sufficient for 
such studies. Experimentally soliton like initial states 
can be created and their dynamics can be studied us- 
ing single site resolution and addressability available in 
current experiments with ultracold atoms in optical lat- 
tices. Density modulation can be achieved by letting the 
system reach equilibrium with a specified inhomogeneous 
potential. Phase modulation can be done by applying a 
short potential pulse. We note that clear signatures of 
soliton dynamics can be achieved by starting with ini- 
tial states that have density modulation and no phase 



FIG. 7: Soliton dynamics in the presence of parabolic poten- 
tial. Soliton propagation is stable when parabolic potential 
is not taking the system across boundaries of different re- 
gions(figs a-c). Crossing the phase boundary causes the soli- 
ton break-up (fig d). Green line marks the boundary between 
regions 1+ and 11+ . 



imprinting. In this case a pair of solitons appears and 
propagates in the opposite directions. This process takes 
place on relatively short time scales accessible to exper- 
iments. Additional details about this procedure are dis- 
cussed in the supplementary material. We also note that 
to observe soliton dynamics one does not need to cre- 
ate initial configuration, which match theoretically cal- 
culated solitonic solutions exactly Numerical analysis 
of (2) shows that initial conditions with spatially local- 
ized excess or deficit of magnetization (density) separate 
reasonably fast into solitons and the wave background. 
On the other hand observing all stages of the magnetiza- 
tion (density) step decay requires longer times and larger 
system sizes. Another important experimental consider- 
ation, which we did not address so far, is the presence of 
the parabolic confining potential in all of the currently 
available experimental set-ups. In fig 7 we show that as 
long as the chemical potential change is not taking the 
system across one of the boundaries in fig. 1, it has no 
strong effect on the soliton propagation. However cross- 
ing any of the boundaries leads to the soliton break-up. 

Summary. We analyzed semiclassical dynamics of 
anisotropic Hciscnbcrg models in dimension higher than 
one. Combining analytical study of small amplitude, 
long wavelength excitations with numerical studies of 
large amplitude, short wavelength excitations we demon- 
strated the existence of a universal dynamical phase di- 
agram, in which different regions can be distinguished 
based on the character of one dimensional solitonic exci- 
tations, and stability of such solitons to transverse mod- 
ulation. Universality of dynamics, which we find from 
direct solution of lattice equations, is very intriguing. 
Our model is not a special lattice regularization of an 
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intcgrablc continuum system. Wc analyze lattice model 
describing real physical systems which, in principle, can 
contain many terms breaking universality. 
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Supplementary material 

Semiclas steal dynamics. Wc consider time-dependent 
variational wavefunctions 



i*w> = n 



0i(t) 6i(t) mm 

sin e ^ | l) i + cos e * | f), 



(3) 

To project quantum dynamics into these wavefunctions 
we define the Lagrangian [31, 50, 51] 

C = -i(*|4|*) + <*|F|¥) = 



Linearizing equations of motion around a uniform state 
we have 



1.2 



: ±\{r 1 ,r 2 )r 



1.2 



(7) 

with A = 4 v /2J L (l-^i 2 )( 71 - J z ). Here r 1 ' 2 are 
Riemann invariants [31] that have physical interpretation 
of left and right moving perturbations. When A is real we 
have hyperbolic regime of the system, when A is imag- 
inary we have elliptic regime. Special lines J 2 /Jj_ = 1 
and fig = ±1 have A = and correspond to the parabolic 
regime. 

Relation to Kortewegde Vries equations and solitons. 
To understand the character of solitons in the hyperbolic 
regime we take equations of motion obtained from (6) 
and keep only the lowest order terms in nonlinearity. We 
further assume that left and right moving components of 
the perturbation can be considered separately (e.g. we 
assume that left/right moving parts separate spatially 
before effects of nonlinearity become important) 



= -6lM»/J±(J±-J*) rl r x + 



(8) 



h2 j 2J x (l-t4) 



J_i — J z 



Mo 



Mo 



' XXX 



cost/,- COS( 



(4) 



6 Mo V J±(J± ~ Jz) r 2 r 



x 



(9) 



_ J ± Ei,i=j±l Sm <?i Sin ^ C0S(V? 4 - tfij) 

and write Lagrange equations for all variables. We obtain 

(pi sin 0i = 4 J z sin 0i (cos 9i+i + cos 0j_i) — 
-4Jj_ cosOi (sin(9 i+ i cos(^.j - (fi+i) + sin0i_i cos(^i - <pi 

Oi = -4J_l (sin^+i sin(<y9j - tp l+ i) + sin^-i sin(^ - 

.(5) 

Introducing variables Xi, yi and Zi we obtain equations 
(2) of the main text. 

Nonlinear Hydrodynamics. To describe the regime of 
long wavelength modulations wc introduce slow vari- 
ables in space and time X = hx, T = ht, where h 
is the lattice constant. We also define ^ = cos 9 and 
a(X, T) = hip(X,t). Taking h to be a small parameter 
we obtain 

C = \ a T /U - 2 Jj_ (1 - /i 2 ) cos ax - 2 J z /i 2 + 
+ h 2 Jl (1 - fi 2 ) (\a X xx sinox + \ a 2 xx cosct y ) + 



J±. — Jz 



/'0 



6 1-m 2 



qJz ) r XX 



Equations (8), (9) are of the KdV type except for special 



mes 



l)) 



h 2 J ± ^ cos a x - h 2 J zl ni X x + 0(h 4 



(6) 



Ho = 



,'1 Mo 



- - J z ^ (10) 
6 



Analysis of solitons, including their stability to transverse 
modulation can now be performed as discussed in ref. 
[31]. Equation (10) gives phase boundaries shown in fig. 
1 of the main text. 

Imperfect preparation of the initial state. 

An important question for experimental realizations 
of solitons is dynamics from initial states that do not 
match exactly soliton solutions. Fig. 8 shows that in this 
case the system separates perturbations into solitons and 
wave-trains. 

Parabolic regime. 

In Fig. 9 we show dynamics of solutions in the 
parabolic regime separating the elliptic and hyperbolic 
regions of the phase diagram 1 in the main text. 
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FIG. 8: Dynamics from a general initial state, which does 
not match a solitonic solution. Red lines show magnetization 
(density), blue lines show superfluid velocity, <f>i — <f>i+i. Initial 
parameters correspond to region I+. Different plots are for 
tj± = 0, 5, 10, 15. Note separation of perturbations into 
solitons and wave-trains. 
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